clear all;

n=1000000;
random=rand(n,10);

cd ..; cd ..; cd ..; cd ..; 
load('code_main_model\6_implied_w_r\computation_results\data.mat',...
    'nmkt','pdf_A_forward','gridb','gridb_pdf','gridtheta','gridtheta_pdf','gridz','trans_z','nz','nb','nb_pdf','ntheta','ntheta_pdf','gb_backward','gtheta_backward','gk_backward','gl_backward','adjust_backward','nu','alpha','period_trans');
cd 'code_figure_table\maintext\Table 2\model';
load('rev_backward.mat','rev_backward');

gb=gb_backward;
gtheta=gtheta_backward;
gk=gk_backward;
gl=gl_backward;
grev=rev_backward;

pdf=squeeze(pdf_A_forward(1,:,:,:,:));

% turn pdf to agent
num_people=n*pdf;
num_people=round(num_people);
nagent=sum(num_people,'all');

% for each simulated agent, keep track the market, z, b, theta
mkt(1:nagent,1:11)=0;
b(1:nagent,1:11)=0;
z(1:nagent,1:11)=0;
theta(1:nagent,1:11)=0;
k(1:nagent,1:11)=0;
l(1:nagent,1:11)=0;
rev(1:nagent,1:11)=0;

n_count=0;
for i_mkt=1:nmkt
    for i_z=1:nz
        for i_b=1:nb_pdf
            for i_theta=1:ntheta_pdf
                n_assign=num_people(i_mkt,i_z,i_b,i_theta);
                mkt(n_count+1:n_count+n_assign,1)=i_mkt;
                b(n_count+1:n_count+n_assign,1)=gridb_pdf(i_b);
                z(n_count+1:n_count+n_assign,1)=i_z;
                theta(n_count+1:n_count+n_assign,1)=gridtheta_pdf(i_theta);
                n_count=n_count+n_assign;
            end
        end
    end
end
for i_t=2:11
    mkt(:,i_t)=mkt(:,1); % assume no migration
end

% get talent
cum_trans_z=cumsum(trans_z')';
cum_trans_z(:,end)=1;

random=random(1:nagent,:);
for i_agent=1:nagent
    for i_t=1:10
        z(i_agent,i_t+1)=find(cum_trans_z(z(i_agent,i_t),:)>random(i_agent,i_t),1);
    end
end

% simulation
for i_agent=1:nagent
    for i_t=1:10
        k(i_agent,i_t)=lin_interpo2(nb,ntheta,gridb,gridtheta,squeeze(gk(i_t,mkt(i_agent,i_t),z(i_agent,i_t),:,:)),b(i_agent,i_t),theta(i_agent,i_t));
        l(i_agent,i_t)=lin_interpo2(nb,ntheta,gridb,gridtheta,squeeze(gl(i_t,mkt(i_agent,i_t),z(i_agent,i_t),:,:)),b(i_agent,i_t),theta(i_agent,i_t));
        b(i_agent,i_t+1)=min(gridb(end),lin_interpo2(nb,ntheta,gridb,gridtheta,squeeze(gb(i_t,mkt(i_agent,i_t),z(i_agent,i_t),:,:)),b(i_agent,i_t),theta(i_agent,i_t)));
        theta(i_agent,i_t+1)=min(gridtheta(end),lin_interpo2(nb,ntheta,gridb,gridtheta,squeeze(gtheta(i_t,mkt(i_agent,i_t),z(i_agent,i_t),:,:)),b(i_agent,i_t),theta(i_agent,i_t)));
        rev(i_agent,i_t)=lin_interpo2(nb,ntheta,gridb,gridtheta,squeeze(grev(i_t,mkt(i_agent,i_t),z(i_agent,i_t),:,:)),b(i_agent,i_t),theta(i_agent,i_t));
    end
end

save('rev.mat','nagent','rev')